//
//  LinAlg.h
//  Vector_Math
//
//  Created by John R Peterson on 3/10/12.
//  Copyright 2012 All rights reserved.
//

#ifndef Vector_Math_LinAlg_h
#define Vector_Math_LinAlg_h

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdarg.h>

#define TRUE 1
#define FALSE 0
#define PI 3.14159265
#define BASESIZE 10

typedef struct Vec_struct {
    float *val;
    int size;
    short isColumn;
    short mode; // are the contents being automatically managed, tells free to do the right thing
} Vec;

typedef struct Matrix_struct {
    float **val;
    int rows;
    int columns;
    short mode; // are the contents being automatically managed
} Matrix;

typedef struct VecList_struct {
    int size;
    int maxsize;
    Vec ** list;
} VecList;

/*
 * March 13
 * 
 * It is important to note that all opperations unless explicitly denoted with an I appended to the command
 * will return pointers to new matrix and vector structs, which must be destroyed using Vec_destroy or Matrix_destroy
 * as appropriate to avoid memory leaks
 *
 * Note that I at then end of an operation means inplace, as in one of the arguments is modified
 * these inplace operations return the target for convenience
 *
 * Other functions do not modify input arguments and instead return a new struct
 *
 * Note that there are 4 different versions of multiplication.  
 * Vec_mulitply for two vectors
 * Matrix_muliplyAB for matricies
 * Matrix_multiplyAV for a Matrix times a vector
 * Matrix_multiplyVB for a vector times a matrix
 * To facilitate chaining operations together, multiplication functions that destroy their arguments have been added
 * Matrix_multiplyAB leaves both of its arguments intact, Matrix_multiplyABD computes the answer and then destroys its first argument
 * Matrix_multiplyAB_d destroys its second argument, and Matrix_multiplyABDD destroys both arguments
 * The other multiplication functions follow a similar naming scheme
 *
 * Most functions if given incompatible vectors ie dimensions don't match up will exit the program
 * equals is the exception and will merely return false
 *
 * To use existing arrays without copying the array with Vector math simply use Vec_createMan, but be careful
 *
 * The R appended to the vector_list functions indicates that they will destroy the vec that you are removing from the list
 */

/******************************************************/
// Vector Function Prototypes 
Vec* Vec_add(Vec *a, Vec *b);
Vec* Vec_addI(Vec *target, Vec *b);
Vec* Vec_addID(Vec *target, Vec *b);
Vec* Vec_addD(Vec *a, Vec *b);
Vec* Vec_add_D(Vec *a, Vec *b);
Vec* Vec_addDD(Vec *a, Vec *b);

Vec* Vec_convolveC(Vec* f, Vec *k);             // convolution of discrete function in 1 D crops to length of 1st argument
Vec* Vec_maxConvolveC(Vec *v, Vec *k);          // max convolution useful for dialation, must be used with the appropriate kernel

Vec* Vec_copy(Vec *v);                          // returns a copy of the given vector

Vec* Vec_create(float values[], int size);      // copies the values and returns a new vector
Vec* Vec_createE(int size);                     // returns a new vector with allocated space but does not initialize values
Vec* Vec_createGaussian(float sigma, int size); // create a discrete valued gaussian centered in the middle of the vector
Vec* Vec_createMan(float values[], int size);   // does not copy values, only uses a pointer to it, use with caution
Vec* Vec_createZeros(int size);                 // initializes a vector of zeros

Vec* Vec_crossProduct(Vec *a, Vec *b);
Vec* Vec_crossProductD(Vec *a, Vec*b);
Vec* Vec_crossProduct_D(Vec *a, Vec*b);
Vec* Vec_crossProductDD(Vec *a, Vec*b);

void Vec_destroy(Vec *target);                  // frees associated memory

float Vec_dotProduct(Vec *a, Vec *b);
float Vec_dotProductD(Vec *a, Vec *b);
float Vec_dotProduct_D(Vec *a, Vec *b);
float Vec_dotProductDD(Vec *a, Vec *b);

short Vec_equals(Vec *a, Vec *b);               // are two vectors equal?

float Vec_L2norm(Vec *v);
float Vec_L2normD(Vec *v);                      // destroys input vector
float Vec_L2norm2(Vec *v);                      // squared L2 norm

Matrix * Vec_multiply(Vec *a, Vec *b);
Matrix* Vec_multiplyD(Vec *a, Vec *b);          // after computing the answer destroys the first input
Matrix* Vec_multiply_D(Vec *a, Vec *b);         // after computing the answer destroys the second input
Matrix* Vec_multiplyDD(Vec *a, Vec *b);         // after computing the answer destroys both inputs

Vec* Vec_normalize(Vec *v);                     // returns a normalized vector such that its L2 norm is 1
Vec* Vec_normalizeI(Vec *target);

void Vec_print(Vec *target);

Vec* Vec_scale(Vec *v, float scaleFactor);
Vec* Vec_scaleI(Vec *target, float scaleFactor);

Vec* Vec_subtract(Vec *a, Vec *b);
Vec* Vec_subtractI(Vec *target, Vec *b);
Vec* Vec_subtractID(Vec *target, Vec *b);
Vec* Vec_subtractD(Vec *a, Vec *b);
Vec* Vec_subtract_D(Vec *a, Vec *b);
Vec* Vec_subtractDD(Vec *a, Vec *b);

Matrix* Vec_toMatrix(Vec *v);

Vec* Vec_2Dto3DI(Vec *target);                       // convert from 2D to 3D
Vec* Vec_3Dto2DI(Vec *target);                       // convert from 3D to 2D

Vec* Vec_toHI(Vec *target);                     // converts a vector to homogeneous coordinates inplace
Vec* Vec_RHI(Vec *target);                      // converts a vector from homogeneous coordinates to normal in place

Vec* Vec_transpose(Vec* v);                     // returns a new vector 
Vec* Vec_transposeI(Vec *target);               // changes parameter in existing vector

/******************************************************/
// Matrix Function Prototypes
Matrix* Matrix_add(Matrix *A, Matrix *B);
Matrix* Matrix_addI(Matrix *target, Matrix *B);
Matrix* Matrix_addID(Matrix *target, Matrix *B);
Matrix* Matrix_addD(Matrix *A, Matrix *B);
Matrix* Matrix_add_D(Matrix *A, Matrix *B);
Matrix* Matrix_addDD(Matrix *A, Matrix *B);

Matrix* Matrix_convolveC(Matrix *F, Matrix *K);     // 2D convolution, crops to size of first argument, 2nd argument is kernel, must be odd
Matrix* Matrix_maxConvolveC(Matrix *F, Matrix *K);  // 2D max convolution, must be provided an appropriate kernel

Matrix* Matrix_copy(Matrix *m);

Matrix* Matrix_create(int rows, int columns, ...);              // expects # of rows comma seperated float *s of length columns
Matrix* Matrix_createCircle(float radius);                      // creates a circle kernel for use with maxconvolution
Matrix* Matrix_createDiag(float *diag, int rows, int columns);  // creates matrix with diag on the diagonal of the matrix
Matrix* Matrix_createV(float *values, int rows, int columns);   // creates a matrix from a single vector of length rows * columns, assumed to be layed out [row1, row2, ...]
Matrix* Matrix_createE(int rows, int columns);                  // creates uninitialized matrix, pointers are all set up
Matrix* Matrix_createGaussian(float sigma, int size);           // creates a 2D discrete gaussian kernel
Matrix* Matrix_createIdentity(int rows, int columns);           // creates identity matrix
Matrix* Matrix_createMan(float **values, int rows, int columns);// good luck, you have to set up your pointers yourself and everything
Matrix* Matrix_createRBT2D(float dx, float dy, float dtheta);   // rigid body transformation matrix in 2D
Matrix* Matrix_createZeros(int rows, int columns);              // initializes all elements to zero

Matrix* Matrix_rotateX(float theta);                            // 3D homogeneous rotation matricies
Matrix* Matrix_rotateY(float theta);
Matrix* Matrix_rotateZ(float theta);
Matrix* Matrix_translateX(float x);
Matrix* Matrix_translateY(float y);
Matrix* Matrix_translateZ(float z);
Matrix* Matrix_translate(float x, float y, float z);

void Matrix_destroy(Matrix *m);                                 // free the struct

float Matrix_det(Matrix *m);                                    // matrix determinant 2x2 or 3x3 only

void Matrix_eigDecomp(Matrix *A, Matrix **ans);                 // performs eigenvalue decomposition 2x2 only places resulting matricies in ans first return is a diagonal matrix of eigenvalues, 2nd is a matrix whose columns are the associated eignevectors

short Matrix_equals(Matrix *A, Matrix *B);        

Vec* Matrix_getRow(Matrix *m,int row);
Vec* Matrix_getColumn(Matrix *m, int column);

Matrix* Matrix_inverse(Matrix *m);                      // performs matrix inversion for 2x2 or 3x3 case only

Matrix* Matrix_multiplyAB(Matrix *A, Matrix *B);        // performs matrix multiplication with another matrix
Matrix* Matrix_multiplyABD(Matrix *A, Matrix *B);       // performs multiplication and destroys 1st argument
Matrix* Matrix_multiplyAB_D(Matrix *A, Matrix *B);      // performs multiplication and destroys 2nd argument
Matrix* Matrix_multiplyABDD(Matrix *A, Matrix *B);      // performs multiplication and destroys both arguments

Vec* Matrix_multiplyAV(Matrix *A, Vec *v);              // with a vector on the right and outputs a vector
Vec* Matrix_multiplyAVD(Matrix *A, Vec *v);             // destroys 1st argument
Vec* Matrix_multiplyAV_D(Matrix *A, Vec *v);            // destroys 2nd argument
Vec* Matrix_multiplyAVDD(Matrix *A, Vec *v);            // destroys both arguments

Vec* Matrix_multiplyVB(Vec *v, Matrix *B);              // with a vector on the left
Vec* Matrix_multiplyVBD(Vec *v, Matrix *B);             // destroys 1st argument
Vec* Matrix_multiplyVB_D(Vec *v, Matrix *B);            // destroys 2nd argument
Vec* Matrix_multiplyVBDD(Vec *v, Matrix *B);            // destroys both arguments

void Matrix_print(Matrix *m);

Matrix* Matrix_scale(Matrix *m, float scaleFactor);
Matrix* Matrix_scaleI(Matrix *target, float scaleFactor);

Matrix* Matrix_subtract(Matrix *A, Matrix *B);
Matrix* Matrix_subtractI(Matrix *target, Matrix *B);
Matrix* Matrix_subtractID(Matrix *target, Matrix *B);
Matrix* Matrix_subtractD(Matrix *A, Matrix *B);
Matrix* Matrix_subtract_D(Matrix *A, Matrix *B);
Matrix* Matrix_subtractDD(Matrix *A, Matrix *B);

Matrix* Matrix_transpose(Matrix *m);

/******************************************************/
// Vector List function prototypes - like array lists but not as good
void VecList_add(VecList *vlist, Vec *v);             // adds the selected vector from the list
void VecList_combineI(VecList **targetprt, VecList **blistprt); // appends blist to target and destroys blist
VecList* VecList_copy(VecList *that);
VecList* VecList_create(void);
VecList* VecList_createS(int initSize);               // creates an empty list with the given initial size
void VecList_destroy(VecList *vlist);                 // does not destroy the vectors contained within
void VecList_destroyR(VecList *vlist);                // destroys each vector it points to
Vec* V_getV(VecList *vlist, int index);                 // returns pointer to vector
short VecList_remove(VecList *vlist, int index);      // removes the vector from the selected index returns true if successful, does not destroy the removed vector
short VecList_removeR(VecList *vlist, int index);

#endif
